#include "Result.h"
//calculate the results wanted**************************************************************************
bool RCS::GetResults()
{	
	real_data Theta_Step = ((1 == Theta_Num) ? 0.0:(Theta_End-Theta_Start)/(Theta_Num-1));
	real_data Phi_Step = ((1 == Phi_Num) ? 0.0:(Phi_End-Phi_Start)/(Phi_Num-1));
	blitz::Array<real_data,2> FarField(Theta_Num,Phi_Num);
	real_data max = 0;
	for (int theta=0; theta<Theta_Num; theta++)
	{
		for (int phi=0; phi<Phi_Num; phi++)
		{
			FarField(theta,phi) = Farfld((Theta_Start+theta*Theta_Step)/180.0*PI,(Phi_Start+phi*Phi_Step)/180.0*PI,Polar_Type);
			if (FarField(theta,phi)>max) max = FarField(theta,phi);
		}
	}
	std::ofstream fout("RCS.dat");
	for (int theta=0; theta<Theta_Num; theta++)
	{
		for (int phi=0; phi<Phi_Num; phi++)
		{
			fout.setf(std::ios_base::scientific,std::ios_base::floatfield);
			fout.setf(std::ios_base::showpoint);
			fout.setf(std::ios_base::right,std::ios_base::adjustfield);
			real_data Mag = 4.0 * PI * FarField(theta,phi)* FarField(theta,phi);
			fout<<10.0 * log10(Mag+1e-30)<<std::endl;
			// 			fout<<FarField(theta,phi)<<endl;
		}
	}
	fout.close();
	return true;
}

real_data RCS::Farfld(real_data theta, real_data phi, int PolarType)
{
	real_data thesin = sin(theta);
	real_data thecos = cos(theta);
	real_data phisin = sin(phi);
	real_data phicos = cos(phi);
	Position3D Pol, Rk;
	Rk = Position3D(Ex->Wave_Num * thesin * phicos, Ex->Wave_Num * thesin * phisin, Ex->Wave_Num * thecos);
	if (Polar_Type == 0)
	{
		Pol = Position3D(thecos * phicos, thecos * phisin, -thesin);
	}
	else if (Polar_Type == 1)
	{
		Pol = Position3D(-phisin, phicos, 0.0);
	}
	Complex3D MoM_Far(0.0,0.0,0.0);
	for (int i=0; i<Geo->Edge_Num; i++)
	{
		int node1 = Geo->Edge(i)[0];
		int node2 = Geo->Edge(i)[1];
		Complex3D tmp(0.0,0.0,0.0);
		for (int LL=0; LL<2; LL++)
		{
			int pat = Geo->Edge(i)[LL+2];
			int node3 = IGetN3(pat,node1,node2);
			real_data signal(1.0-2*LL);

			Position3D Rho(Geo->Patch_Ctr(pat)-Geo->Node(node3));
			tmp +=Rho * signal * exp(dComplex(0.0,1.0)*dot(Rk,Geo->Patch_Ctr(pat)));
		}
		MoM_Far +=  Geo->Edge_Length(i) * Geo->I(i) * tmp;
	}
	MoM_Far = MoM_Far * 0.5;
	MoM_Far = -dComplex(0.0,1.0) * Ex->Omega * Mu/(4*PI) * (MoM_Far);
	return abs(dot(MoM_Far,Pol));
}
